Elasticity of Stiff Polymer Networks 
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We study the elasticity of a two-dimensional random network of rigid rods ("Mikado model"). 
The essential features incorporated into the model are the anisotropic elasticity of the rods and the 
random geometry of the network. We show that there are three distinct scaling regimes, character- 
ized by two distinct length scales on the elastic backbone. In addition to a critical rigidiy percolation 
region and a homogeneously elastic regime we find a novel intermediate scaling regime, where the 
elasticity is dominated by bending deformations. 

PACS numbers: 87.16.Ka, 62.20.Dc, 82.35.Pq 



The elasticity of cells is governed by the cytoskele- 
ton, a partially crosslinked network of relatively stiff fil- 
aments forming a several 100 nm thick shell called the 
actin cortex pj . While the statistical properties of single 
cytoskelctal filaments arc by now relatively well under- 
stood 0, Q , theoretical concepts for the elasticity of stiff 
polymer networks are still evolving. One major open 
question is to understand how stresses and strains are 
transmitted in such networks. In synthetic gels that are 
formed by rather flexible chain molecules the response 
to macroscopic external forces is - on the level of single 
filaments - isotropic and entropic in origin. It is gen- 
erally believed that macroscopic stresses are transmit- 
ted in such a way that the local deformations within the 
network stay affine, i.e. that the end-to-end distance of 
individual filaments follows the macroscopic shear defor- 
mation 0| . In contrast, the building blocks of the actin 
cortex are semiflexible polymers, whose hallmark is an 
extremly long persistence length £ p , which is comparable 
to the total contour length £. As a consequence, the re- 
sponse of such stiff polymers to external forces shows a 
pronounced anisotropy [5j- Consider a semiflexible poly- 
mer with one end clamped at a fixed orientation. When 
forces are applied at the other end transverse to the tan- 
gent vector at the clamped end, the response may be 
characterized by a transverse spring coefficient k±(£) — 
3k/£ 3 proportional to the bending modulus k. Whereas 
this response is of purely mechanical origin, the linear 
response for longitudinal forces is due to the presence 
of thermal undulations, which tilt parts of the polymer 
contour with respect to the force direction. The corre- 
sponding effective spring coefficient k\\ (£) = 6k 2 / (k,BT£ A ) 
is proportional to k 2 /T indicating the breakdown of lin- 
ear response for very stiff filaments. In a typical network 
one expects the distance between crosslinks £ c to be much 
smaller than the persistence length and filament length. 
Hence we have k\\(£ c )/k±(£ c ) = 2£ p /£ c 3> 1, i.e. the elas- 
tic response of the filaments is indeed highly anisotropic. 

These anisotropic elastic properties of individual fila- 
ments suggests that the macroscopic elasticty of networks 
will not only depend on the number of crosslinks and the 



density of filaments, but also on the geometry and archi- 
tecture of the network. For very regular networks such as 
a triangular lattice the longitudinal spring coefficient k\\ 
will dominate the macroscopic moduli since the net- 
work can not be strained without a change of the end- 
to-end distance of individual polymers. In other regular 
network architectures, the softer bending modes would 
be dominant [jj. Naturally, this will lead to a very dif- 
ferent prediction for the elastic modulus of the network. 
It is not at all obvious what type of network geometry 
(elongation dominated versus bending dominated) is rel- 
evant in less ideal structures with a significant amount of 
disorder such as in typical cytoskeletal networks. 

As a first step towards understanding the elasticity 
of stiff polymer networks we consider a two-dimensional 
model defined as follows (see Fig. We generate the 




FIG. 1: Typical networks at low and high density. Dangling 
bonds, not contributing to the elasticity, have been cut off. 
The stress distribution is shown in false colors; the load on 
a filament increases from blue to red. The left picture is for 
p — 10, system size L = 10, and an aspect ratio a = 0.0001. 
99.99% of the strain energy stored in bending modes. In con- 
trast, the right picture shows a network for p — 50, L = 2, and 
a = 0.01, where only 5% of the strain energy is in bending 
modes; the remainder is stored in compression modes. For 
the choice of units see the main text. 

random network by placing N line-like objects of equal 
length I on a plane with area A = L 2 such that both 
position and orientation of the filaments are uniformly 
randomly distributed. Periodic boundary conditions in 
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both directions are used. Upon increasing the line den- 
sity p = N£/A there is a critical threshold p c for geo- 
metric percolation Numerical simulations 9j show 
that the correlation length £ ~ (p — p c )~ u of the in- 
cipient infinite percolation cluster scales with a critical 
exponent v = 4/3, identical to the value obtained for 
random site percolation on a lattice |ld| . Transport of 
scalar quantities like the conductivity is also in the same 
universality class as lattice models . In order to study 
the transport of non-scalar quantities like shear stress we 
need to specify how forces are transmitted between the 
building blocks of the network. In our "Mikado model" 
the building blocks are homogeneous elastic rods char- 
acterized by a Young modulus E and a circular cross- 
section of radius r. Wherever two sticks intersect they 
are connected by a crosslink with zero extensibility. In 
the cytoskeleton one finds a variety of linker proteins 
with a range of mechanical properties . Here we re- 
strict ourselves to crosslinks that either fix the relative 
orientation of the rods ("stiff crosslinks") or allow free 
rotation ("free hinges"). Similar to thermally fluctuat- 
ing semiflexible polymers, the elastic response of a stick 
segment between two neighboring crosslinks is character- 
ized by length dependent force constants for compres- 
sion or elongation, k comp (£ c ) = irr 2 E/£ c , and bending, 
fc b cnd(4) = k±(£ c ) = {3/4)irr A E/£ 3 c . The distance be- 
tween two crosslinks £ c shows a Poissonian distribution, 
where the average distance of crosslinks along a filament 
scales as the inverse of the line density, £ c = ir/p ||. 
While this is a purely mechanical model, it still captures 
the essential feature that for typical densities of the net- 
work the compressional stiffness is much larger than the 
bending stiffness, k comp (£ c ) / k hcnd (£ c ) = (4/3)^/ r 2 > 1. 
It leaves out steric effects due to thermal fluctuations of 
the filaments, which give rise to the plateau modulus in 
solutions 0. 

Consider the energy of the network as a function of the 
deviations of the positions of all intersection points and 
the rod orientations at the intersection points from their 
initial values. For small deformations of the network, this 
function can be approximated by a quadratic form that 
vanishes for vanishing deviations, as — by construction — 
the undeformed network is not prestressed. To analyse 
the elastic properties of the model network, a shear de- 
formation respecting the periodic boundary conditions 
is enforced by demanding that corresponding points on 
the left and right boundary of the simulation cell un- 
dergo equal displacements while the displacements of cor- 
responding points on the upper and lower boundary of 
the cell must agree vertically but differ horizontally by a 
distance A = jL, where 7 is the shear strain. The orien- 
tation of the rods at corresponding points on the bound- 
ary are required to be equal. The remaining degrees of 
freedom are then allowed to relax, i.e. the harmonic ap- 
proximation to the energy of the network is minimized 
in the presence of the constraints. The derivative of the 



resulting energy of the deformed state with respect to the 
strain 7 is proportional to the shear modulus. In pinciple, 
this reduces the determination of the shear modulus of a 
given network to the solution of a linear equation. How- 
ever, for interesting parameters (thin rods), the problem 
is numerically highly unstable as we are searching for 
the lowest point of a complicated high-dimensional val- 
ley with extremely steep slopes but hardly varying base 
altitude. This problem is best left to one of the com- 
mercially available finite element solvers which have seen 
many years of careful optimization and testing. The re- 
sults presented below were obtained using the program 
Nastran by MSC Software. 

In the following discussion we take the rod length £ 
as unit length and k/£ 3 as unit for the elastic modulus. 
Then the independent parameters are the densitiy p, the 
system size L and the aspect ratio a — r/£ of the rods. 
Note that the latter is a measure of the relative magni- 
tude of compressional to bending stiffness. 

We start with an analysis of the elasticity close to the 
percolation threshold. For stiff crosslinks we find that 
the percolation threshold is the same for rigidity as for 
connectivity percolation, p c — 5.71. For free hinges a 
higher line density of filaments is needed, p c — 6.7, for 
the network to become rigid. This agrees well with re- 
cent results, p c — 6.68, for stiff fiber networks 0], where 
the crosslinks are fixed in space but the angles between 
the fibers can vary. In both cases, we find that the 
shear modulus G vanishes as the line density of sticks 
approaches the critical value p c , according to a power 
law G ~ (p — Pc)' 1 ■ For our numerical analysis on a finite 
lattice we expect the shear modulus to obey the following 
finite size scaling law 

G = L-^h(L/0, (1) 

where the scaling function behaves as h{x) ~ x^' 11 and 
h(x) ~ 1 for large and small values of the scaling vari- 
able x = L/£, respectively. Fig. |21 shows that the data 
collapse works very well for densities ranging from values 
close to the percolation threshold p c up to p « 20. For 
the data shown, L ranges from 2 to 30. For larger den- 
sities, systematic deviations are clearly visible. This will 
turn out to be a very interesting observation, as we will 
discuss in detail below. We get the best data collapse 
in the critical region if we choose the values 2.4 ± 0.2 
and 2.3 ± 0.2 for the critical exponent p,jv in the case 
of stiff crosslinks and free hinges, respectively. Since the 
difference between the exponents is within the statisti- 
cal error, we can make no definite conclusion whether 
networks with free hinges and stiff crosslinks belong to 
different universality classes for elasticity percolation. 

The rigidity exponent p w 3.15 ± 0.2 is significantly 
lower than in other classes of continuum percolation mod- 
els, like the "Swiss-cheese model", where /i«5 [lfil Hfi| . 
It is also lower than the value u « 4 for lattice models 
with bond-bending forces jlOl Il7j| . Hence it seems likely 
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FIG. 2: Double logarithmic plot of the scaling function h(x) 
for the shear modulus of the "Mikado model" with stiff 
crosslinks (top) and free hinges (bottom) as a function of 
x — L\5p\ u with 8p — p — p c for a series of densities p in- 
dicated in the graphs. Note that for finite systems the shear 
modulus is also nonzero below p c (lower branches in both 
plots) . 



that the "Mikado model" constitutes a new universality 
class for rigidity percolation. Similar results have inde- 
pendently been found in Ref. [is| . 

Now we come back to the above mentioned system- 
atic deviations from the scaling law, Eq. ^ a t densities 
above p w 20. To understand these better, let us have 
a closer look at the shear modulus as a function of the 
elastic moduli of individual filaments for densities not 
too close to the percolation threshold. In this regime 
the shear modulus becomes independent of system size 
for moderately large systems; for the following results 
we have chosen systems satisfying L/£ > 200. Fig. |3| 
shows the shear modulus as a function of a for a se- 
ries of densities; we have communicated a preliminary 
version of these data in Ref. 0. Note that febendW 
is effectively kept constant since we are measuring all 



10 15 



10 10 



>■-■ 

V 

•3- 




10 



FIG. 3: Double logartithmic plot of the shear modulus G as 
a function of a for fixed fcbend(^)- Data are shown for free 
hinges. 

ingly different regimes. For high densities and/or thick 
rods (a > 0.1), where compressional stiffness is lower or 
comparable to the bending stiffness (lower right part of 
Fig. I3J , the shear modulus scales linearly with the fila- 
ment compressional modulus and the number of filaments 
per unit area, G ~ (p — p c )a~ 2 . Such a linear regime has 
also been found in a series of studies on random fiber 
networks [2(J . It is by now well established that the elas- 
tic modulus can be described quantitatively in terms of 
effective medium models |2l|. Hence, in the high line 
density regime the network behaves as a homogeneously 
elastic medium, dominated by the compressional mod- 
ulus of the individual filaments. As a consequence, lo- 
cal deformations follow a macroscopic shear in an affine 
way. This has to be contrasted with the elastic behav- 
ior for slender rods with low aspect ratios (a m 10~ 5 
for the higher densities), where bending becomes the 
softer mode. Then, one finds an extended plateau re- 
gion, which broadens significantly with lowering the line 
density, where the shear modulus becomes completely in- 
dependent of k comp (£) ~ a~ 2 k hcnd (£) 19]. This strongly 
suggests that in this regime the macroscopic elasticity of 
the network is dominated by bending stiffness of indi- 
vidual filament. This conclusion is corroborated by the 
observation that almost all of the energy stored in the 
deformed network is accounted for by transverse defor- 
mation of the rods (compare Fig.^l. Another remarkable 
feature of this plateau regime is the strong dependence of 
the shear modulus on line density. We find G ~ (p~p c Y 
with a rather large exponent p! « 6.7. From the above 
analysis it may seem as if the anomalous elasticity in 
the plateau regime and the homogeneous elasticity in the 
affine regime are two separate phenomena, and one might 
wonder how one emerges from the other. To analyze this 
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relation we try a crossover scaling ansatz, 

G = ( P - PcY'gWp - Pcf] = r»' /v '~g{<*IO , (2) 

where we have defined a new length scale £' ~ (p— p c )~ v ■ 
For this ansatz to reduce to the modulus expected in the 
affine region, the scaling function g{x) needs to scale as 
g(x) ~ x~ 2 for x 3> 1 and the exponents need to obey 
the scaling relation p! = 2v' + 1. In the plateau regime, 
g(x) is expected to be constant. As shown in Fig. 0] we 
obtain an excellent scaling collapse for over almost eight 
orders of magnitude in the scaling variable x — a/£' us- 
ing v 1 = 2.83 or equivalently p' = 6.67 and the critical 
line density p c w 5.71, associated with connectivity per- 
colation. Additionally, the scaling function g(x) displays 
the expected behavior. Meeting both of these require- 
ments is highly nontrivial, and gives strong evidence for 
the anomalous scaling law in Eq. 
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FIG. 4: Scaling plot of the shear modulus for free hinges for a 
series of densities above p = 15 indicated in the graph (same 
data as in Fig.0. Data collapse to the crossover scaling form, 
Eq. is obtained with v' = 2.83. Note that here and in all 
other figures the unit of length is £ and the unit of the shear 
modulus G is k/£ 3 . 

The existence of such a broad scaling regime far from 
the percolation threshold is a surprising and intriguing 
feature of stiff polymer networks. At the moment we 
are lacking a complete understanding of its physical ori- 
gin. In particular, the geometrical significance of the 
new length scale £' remains unclear. One may speculate 
that the anomalous scaling behavior is a subtle conse- 
quence of the interplay between quenched randomness of 
the network structure and long-range correlation effects 
induced by the stiffness of the filaments. An immediate 
consequence of the scaling form, Eq. |2 is the existence of 
a crossover line density p C ross scaling as £p cm ss ~ or 1 /" , 
where we have re-introduced units of length £. This im- 
plies that increasing filament length at constant line den- 
sity drives the system towards the affine regime, in accord 
with Ref. 185. 



While these results for an idealized two-dimensional 
model are certainly not straightforwardly applicable to 
three-dimensional cytoskeletal networks, one may still try 
to get an idea for the scales involved. We expect that 
network densities can be compared roughly by using the 
average distance £ c between intersections as a measure: A 
cytoskeletal network might have £ c w 0.1 pro. with typical 
filament lengths of 2 pm. These values correspond to a 
two-dimensional line density of p 20 and an aspect 
ratio of a R* 0.002, which would place a typical actin 
network in the bending dominated intermediate regime. 

Understanding the full complexity of cytoskeletal net- 
works certainly merits further theoretical and experimen- 
tal work. Building on the knowledge gained from our 
idealized model, future investigations may among many 
other questions want to address three-dimensional sys- 
tems, polydispersity, thermal fluctuations or even the ki- 
netics of the crosslinking molecules. 

We acknowledge M. Alava and K. Kroy for useful dis- 
cussions and comments, and P. Benetatos for a critical 
reading of the manuscript. 
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